#include "tbparameters.hpp"


/* Optical transition matrix:
    
     We consider weak optical field and construct
     
     U(t) = 1 - i Integral( A(t').V, t'=0...t)

     where A(t) is the vector potential for the optical field and since
     we negelect, at least for now, the dynamics within the (~100 fs) pulse
     we will simply set, t > the time for the pulse to vanish. 

     We then re-construct U to obtain a matrix that is exactly unitary
     
     U = exp( U - 1)

    */
    

void opticalinjection_1storder(double pulsewidth,
		      double centerfreq,
		      cx_vec pol,
		      vec E,
		      cx_mat VelocityMatrix[],
		      cx_mat U){

  int n = E.n_elem;
  cx_mat Y(n,n);
  cx_double matelem;
  double A;

  for (int i=0; i<n; i++){
    for (int j=0; j<n; j++){
      
      matelem = 0.0;
      for (int mu=0; mu<3; mu++)
	matelem += VelocityMatrix[mu](i,j)*pol(mu);

      A = (E[i]-E[j] - centerfreq)*pulsewidth/hbar;
      Y(i,j) = matelem * exp(-A*A*2);
    }
  }

  cx_mat W;
  vec d;
  eig_sym(d,W,Y);
  cx_mat D(n,n);
  D.zeros();
  for (int i=0;i<n;i++)
    D(i,i) = exp(-cx_double(0.0,-1.0)*d(i));

  U  = W*D*W.t();

  return;
}

void opticalUnitaryMatrix(double pulsewidth,
			  double centerfreq,
			  cx_vec pol,
			  vec E,
			  cx_mat VelocityMatrix[],
			  cx_mat U){

  for (int i=0; i<n; i++){
    for (int j=0; j<n; j++){
      
      matelem = 0.0;
      for (int mu=0; mu<3; mu++)
	matelem += VelocityMatrix[mu](i,j)*pol(mu);
      
      A = (E[i]-E[j] - centerfreq)*pulsewidth/hbar;
      Y(i,j) = matelem * exp(-A*A*2);
    }
  }
  
}
